##########################################################################################

library(data.table)
library(optparse)
library(ArchR)
library(Seurat)
library(ggsci)
library(dplyr)
library(GenomicFeatures)

##########################################################################################
option_list <- list(
    make_option(c("--comine_data_file"), type = "character"),
    make_option(c("--nclust"), type = "character"),
    make_option(c("--gene_type"), type = "character"),
    make_option(c("--scriptPath"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 整合atac和rna的文件
    comine_data_file <- "~/20231121_singleMuti/results/qc_atac_v3/germ/testis_combined_peak.combineRNA.qc.Rdata"

    ## 表达文件
    rna_data_file <- "~/20231121_singleMuti/results/qc_atac_v3/germ/testis_combined.annotationCellType.qc.Rdata"

    ## 既往研究整理的代码
    scriptPath <- "~/20231121_singleMuti/scripts/scScalpChromatin"

    ## gtf_file 
    gtf_file <- "~/ref/GTF/20240317_gencode_v32_gene.bed"
    gff3_file <- "~/ref/GTF/gencode.v32.annotation.gff3"

    ## 
    nclust <- 11

    gene_type <- "protein_coding"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/celltype_plot/peak2gene/germ/protein_coding"

}


###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

comine_data_file <- opt$comine_data_file
scriptPath <- opt$scriptPath
nclust <- as.numeric(opt$nclust)
gene_type <- opt$gene_type
out_path <- opt$out_path

dir.create( out_path , recursive = T)

##########################################################################################
## 已发表文献写好的脚本
source(paste0(scriptPath, "/plotting_config.R"))
source(paste0(scriptPath, "/misc_helpers.R"))
source(paste0(scriptPath, "/matrix_helpers.R"))
source(paste0(scriptPath, "/archr_helpers.R"))
source(paste0(scriptPath, "/GO_wrappers.R"))

##########################################################################################
## 导入数据
a <- load(comine_data_file)
# testis_combined_peak_combineRNA
b <- load(rna_data_file)

##########################################################################################

dat_gtf <- fread(gtf_file)
txdb <- makeTxDbFromGFF(gff3_file)

##########################################################################################
## 提取关系的基因集合（编码基因或lcnRNA）
dat_gtf <- subset( dat_gtf , V5 == gene_type )

################################################################################
## 构建基因的位置信息
gene_ensg <- dat_gtf[,c("V6","V4")]
colnames(gene_ensg) <- c("GENEID" , "name")

gene_gr <- genes(txdb, columns = c("GENEID"))

gr_df <- data.frame(
  seqnames = seqnames(gene_gr),
  start = start(gene_gr),
  end = end(gene_gr),
  strand = strand(gene_gr),
  GENEID = names(mcols(gene_gr)$GENEID)
)

## gene name对应多个ensg的，选第一个
gr_df <- merge( gr_df , gene_ensg , by = "GENEID" )
gr_df <- gr_df %>% 
group_by( name ) %>%
summarize( seqnames = seqnames[1] , start = start[1] , end = end[1] , strand = "*"  )

gr <- GRanges(
  seqnames = Rle(gr_df$seqnames),
  ranges = IRanges(start = gr_df$start, end = gr_df$end),
  strand = gr_df$strand,
  name = gr_df$name
)

names(gr) <- gr_df$name

## 保留感兴趣的基因
use_gene <- rownames(scrnat@assays$RNA$counts)[rownames(scrnat@assays$RNA$counts) %in% names(gr)]
gr <- gr[use_gene]

seruat_cds_all <- SummarizedExperiment(scrnat@assays$RNA$counts[use_gene,] , rowRanges = gr)
colnames(seruat_cds_all) <- gsub( "_" , "#" ,  colnames(seruat_cds_all) )

## 替换表达矩阵
projHeme5 <- testis_combined_peak_combineRNA
projHeme5 <- addGeneExpressionMatrix(input = projHeme5, seRNA = seruat_cds_all, force = TRUE)
